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Abstract 

We use Molecular Dynamics combined with Dissipative Particle Dynamics to construct a model of 
a binary mixture where the two species differ only in their dynamic properties (friction coefficients). 
For an asymmetric mixture of slow and fast particles we study the interdiffusion process. The 
relaxation of the composition profile is investigated in terms of its Fourier coefficients. While 
for weak asymmetry we observe Fickian behavior, a strongly asymmetric system exhibits clear 
indications of anomalous diffusion, which occurs in a crossover region between the Cases I (Fickian) 
and II (sharp front moving with constant velocity), and is close to the Case II limit. 

PACS numbers: 66.10.-x,66.10.Cb,02.70.Ns 
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I. INTRODUCTION 



The study of anomalous diffusion phenomena in polymeric materials has been of interest 
for decades. One particular instance of non-Fickian transport is Case II diffusion^, which 
is usually observed in glassy polymers subjected to penetration by a low-molecular weight 
solvent. The most characteristic feature of this phenomenon is the development of a sharp 
front in the concentration profile, which advances linearly in time, and ahead of which the 
penetrant concentration is very low. The dynamics is therefore characterized by a single 
parameter, the velocity of the front. In contrast, Fickian (Case I) diffusion is described 
by the diffusion coefficient and a corresponding scaling of position (and amount of sorbed 
penetrant) with i 1 / 2 . 

Case II diffusion has been intensively investigated experiment ally 2 -^^ and 
theoretically^ i 7 ! 8 ^ 9 ^ ! 11 ^ 2 ! 13 ^ 4 ^ 5 ^ 6 ^ 7 ! 18 , and various microscopic and phenomenological 
models have been proposed to explain the observed features. Among the most popular are 
the models of Thomas and Windle 14 (TW), and of Rossi et a / , 15 i 16 i 17 ; which was recently 
extended by Qian and Taylor—. 

Regardless of the details of the various models, the development of the linearly advancing 
front is rather easy to understand 4 . The two decisive preconditions are (i) the existence of 
a strong disparity in the mobilities of the two pure species (glassy matrix vs. penetrant), 
and (ii) the "plasticizing effect", i. e. a strong enhancement of the mobility of the slow 
species (constituent of the glassy matrix) if its molecules are surrounded by those of the fast 
species (penetrant). As soon as the slow molecules are plasticized at the front, they quickly 
make room for the penetrant, which in turn is rapidly transported into the thus-opened free 
volume. This results in a spatially constant concentration profile behind the front. The rate 
of plasticization (matrix dissolution) is determined by the penetrant concentration at the 
front, and therefore remains constant. This gives rise to a constant- velocity front. 

One should note that this generic scenario does not necessarily require a polymeric matrix 
- the "dynamic contrast" between the two species is not specific to polymers. In fact, 
quite similar behavior has been observed in the dissolution of low molecular weight silicate 
glasses 1 - 9 -, where the constant rate of dissolution was attributed to a constant rate of chemical 
reaction (hydration) at the surface. Nevertheless, most studies are concerned with polymer 
systems, where the phenomenon is of particular practical relevance. There are additional 



2 



polymer-specific effects beyond simple Case II behavior, like entanglements and the resulting 
stresses in the matrb*^, which are technologically important. 

For such reasons, one would like to be able to simulate the phenomenon on the com- 
puter in order to gain insight into the underlying molecular mechanisms. A recent attempt 
has been undertaken by M. Tsige and G. S. Gres t 20 i 21 who investigated a coarse-grained 
polymer-solvent system by means of Molecular Dynamics and grand canonical Monte Carlo 
simulation. However, these authors did not observe any deviation from Fickian behavior, 
probably because the simulations failed to reach the necessary separation of time scales, 
even though the polymers were able to vitrify. 

In the present study we attempt a quite different simulational approach in order to create 
the necessary dynamic disparity. As outlined above, one expects that Case II behavior 
(i. e., a linearly propagating sharp front, and a flat concentration profile behind the front) 
should occur as soon as this disparity is sufficiently strong, and one takes into account the 
"plasticizing" effect - even if the model does not exhibit any typical polymer features. As 
our results indicate (see below), our model, though very simple, seems to grasp the main 
features on non-Fickian diffusion at least qualitatively. We believe that a combination of our 
approach with a polymer model would prove promising in revealing the molecular aspects 
of anomalous diffusion in glassy polymers. The purpose of our study is hence twofold: On 
the one hand, we wish to further corroborate the view that dynamic disparity is indeed the 
decisive ingredient for Case II behavior. On the other hand, we propose a methodological 
advancement in the microscopic simulation of such phenomena. 

The simplest model in this spirit is a binary system of particles which are completely 
identical with respect to their static properties (i. e. their interaction potential), but differ 
in their dynamic properties. One possibility to do this is to run a Molecular Dynamics simu- 
lation of a binary mixture whose species are identical except for their masses. This approach 
has, for instance, been pursued in Ref.— . Another possibility is to dampen the motion of the 
particles by friction and compensate this via stochastic thermal noise (Langevin simulation, 
so-called Stochastic Dynamics). The difference in dynamics would then be implemented via 
different friction constants. One disadvantage of this method, though, is that it does not 
conserve the momentum, and hence does not describe the hydrodynamics correctl y 23 ! 24 . The 
latter might however be important for the phenomenon. Fortunately, there exists a variant 
of Stochastic Dynamics, called Dissipative Particle Dynamics (DPD) 2 ^ 2 ^ , which does not 
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suffer from this deficiency. In this method only relative velocities are dampened so that 
Galilean invariance is restored. Furthermore, the stochastic force acts on pairs of particles 
so that Newton's third law (momentum conservation) holds. For more details, see Sec. |HJ 
We have hence introduced the difference in dynamic properties via suitable DPD friction 
coefficients ( ss , ( s f, and Here, the indices "s" and "f" stand for "slow" and "fast", 
respectively. A friction coefficient ( a! 3 indicates the rate of damping of the relative velocity 
between two particles of species a and (3, respectively. To our knowledge such a "binary gen- 
eralization" of DPD simulations, though fully straightforward, has not yet been considered. 
The advantage of this approach is that such a model is more flexible than just implement- 
ing different masses, since it allows for three parameters rather than just two single-species 
parameters. In particular, it is possible to model the "plasticizing" effect in a very simple 
way by choosing a small value for ( s f (as soon as a slow particle is in a fast environment, its 
motion is no longer dampened). 

In our data analysis, we focus on the composition profiles indicative of the interdiffusion. 
Moreover, instead of just looking at the time scaling behavior of the position of the "in- 
terface" (that is, of the locus of a given composition), we rather compare the profiles (or, 
more precisely, their Fourier components) directly to the solution of the (simple) diffusion 
equation. Fickian and non-Fickian behavior is then detected by determining the scaling of 
various quantities with time. 

The remainder of this article is organized as follows: Section |H] gives a brief outline of 
the standard DPD method. Section IIHI then describes our generalization of the established 
version to the binary case with three friction parameters, and specifies the simulation model 
in detail. Two versions of the model are studied: While Model I exhibits pure Fickian 
diffusion, Model II reveals anomalous Case II behavior. These results are presented in Sec. 
IIVI Finally, Sec. IS concludes with a brief summary. 



II. DISSIPATIVE PARTICLE DYNAMICS 



The equations of motion in the DPD algorithm are 

= —v- —v- = Ft + F ifr) + F { 

at rrii at 
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where r*j is a particle position, p t is a particle momentum, and mi is the mass of particle i. 
Defining = r j — fj = r^r^, (here denotes the respective unit vector), we can write the 
conservative force as 

P l = J2-(dU tJ /dr tJ )r tJ (2) 

3 

that satisfies Newton's third law (here Uij is the interaction potential between particles % 
and j). The friction force F^ r ^ is 

F^^-J^^m-Vj)-^, (3) 

3 

where ( is the friction coefficient. Similarly, we get the stochastic forces along the inter- 
particle axes: 

3 (rt) = £*M (4) 

3 

where a is the noise strength related to the friction coefficient ( through the Fluctuation- 
Dissipation Theorem (FDT) a 2 (r) = ksT£(r) and f]ij{t) is a Gaussian white noise variable 
with the properties 77^ = r] jh (77^) = 0, and (r)ij(t)r) k i(t')) = 2(S ik Sji + 5 u 5 jk )5(t - t'). It is 
then easy to see that the relations 

£ir>=£ir>= (5) 

i i 

hold, i. e. momentum is conserved. 
III. THE MODEL 

The DPD technique allows to construct a simple model of a binary system of slow and fast 
particles. As DPD makes use of friction parameters to control the momentum relaxation and 
the mobility of the particles, it is possible to choose the two types of beads to be completely 
identical as far as their static properties (masses, sizes, and interactions) are concerned. The 
simplifying aspect of such a model is that its equilibrium properties are trivially known - 
the equilibrium state is just a random mixture. 

The equations of motion are similar to those in the standard DPD algorithm. Only the 
friction and the stochastic terms are slightly modified as follows: 

F/i = -Cafiinj) [{Vi -Vj)- fij] fij (6) 
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Fg = v*p(r ij )r Hj (t)r ij (7) 

As functions Ca/3( r ) an d <J a p(r) we simply take constant values up to a certain cutoff for which 
we choose the same value r c as for the interaction potential (see below). The mutual friction 
coefficients ( a/ 3 can take 3 different values: When two slow particles interact ( a/ 3 = ( ss , when 
two fast particle interact ( a/ 3 = Qf, while ( a/ 3 = ( s f when a slow and a fast particle interact. 
cr Q/ g is related to ( a p according to the fluctuation-dissipation theorem: a a p = (2A; B T^ Q/3 ) 1 / 2 . 

The interaction potential is chosen as a repulsive Lennard- Jones potential between 
particles i and j: 

c W -^ l, - ,B - , - , • +i, • ni<r ' m 

[o, Tij>r c 

where r c = 2 1 ^ 6 a LJ is the cutoff radius for all the forces acting on bead i. The unit system is 
defined by setting the particle masses, the energy parameter e, and the length a LJ to unity. 

The initial equilibrium configurations of such a system are prepared in a rectangular box 
(L x x L y x L z , where L x = L y = 36 and L z = 18) with periodic boundary conditions (pbc) 
in x, y and z directions. At time t = a slab spanning half the box (— L x /4 < x < L x /A) 
is defined and the property "slow" is assigned to each particle within that slab. The other 
particles are defined as "fast". Due to the pbc the fast particles (solvent) penetrate the 
region of slow particles (matrix) from both sides and thus create two fronts. 

From this moment on the mixing process is monitored by measuring the density profiles 
of the slow and of the fast particles. The values we choose for the ( aj 3 are: 1) ( ss = 10, 
( s f = 1.0 and Qf = 0.1 (system I), time step 8t — 2 x 10^ 3 for the integration of the 
equations of motion, and 2) ( ss = 1000, ( s f = 1 and £// = 1 (system 11) whereby the large 
friction coefficient ( ss = 1000 requires a very small time step: St = 5 x 10~ 5 . The number of 
particles in the box is 2 x 10 4 at total density p = 0.85 which is kept constant throughout the 
simulation. The starting configurations are fully equilibrated by means of DPD simulations 
using the friction parameter Qf. All the simulation experiments are held at temperature 
UbT = 1.0. We measure the profiles every 10 3 MD steps for system I and every 10 4 for system 
II. 70 independent runs are made for system I up to a time t = 2500 (= 1.25 x 10 6 MD time 
steps) where the final stage of the mixing process is established. The results obtained for 
system II are averaged over 6 independent runs up to time t = 1725 (= 3.45 x 10 7 MD time 
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steps) where the system is still far away from an equilibrium mixture. Note that the time 
step 5t should be carefully chosen since a too large integration step may strongly influence 
the numerics and lead to unphysical results. To this end we have performed additional 
control runs using a five times smaller integration step and verified that our control results 
do not deviate from those which were derived with the reported St. We show a comparison 
of these at the end of the following section. 

IV. SIMULATION RESULTS 

Density profiles obtained by using the adopted simulation method are shown in FigO 
Principally two processes take place during matrix dissolution: 1) the establishment of a 
smeared-out front as a result of Fickian diffusion, and 2) movement of the front position. 
Case II diffusion is observed when the second process is much faster than the first one^. This 
regime is achieved when the difference in the tracer diffusion coefficients of slow particles in 
the matrix D™ and in the solvent D s s is large enough, i. e. D s s 3> D™. This is the assumption 
of Ref.~ while the effect of "plasticizing" in Rossi's models is implemented via D s s = oo. 

The behavior of system I (Fig. [TJi) does not follow this scheme and a deviation from 
Fickian diffusion is not observed. The difference in friction coefficients is small and the 
profile of slow particles in the fast environment looks essentially the same as that of the fast 
particles in the slow matrix. In contrast, a lack of such symmetry, due to the presence of 
both processes, the moulding of the front and its moving in time, can be noticed in system 
II (Fig. m>) where a regime of non-Fickian transport is attained. 

The dynamics of penetration of fast particles in the matrix of slow particles and vice versa 
can be analyzed by looking at the front propagation (the position of each front is determined 
by the inflection point of the respective density profile), Fig. Et, and at the trajectories of 
points with constant density, Fig. |2b- The motion of the front, X in fi, is almost negligible 
in system I (not shown in Fig. |2K)> and only a slight deviation from its constant value at 
early times can be detected. Its velocity is much smaller than the speed of the process of 
front shape formation which governs the behavior of the system. 

A nearly linear dependence of X in fi(t) is observed in system II up to the latest time of 
the simulation, as expected for Case II diffusion. However, also a continuous decay of the 
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front velocity is detected. One could expect that after time t pa 1300 a crossover from Case 
II to Case I behavior develops and the system goes back to Fickian diffusion. 

The three positions of constant density, X i(t), X om (t) and X 0Q1 (t) in system I (Fig. 
I2b) follow the typical behavior of Fickian diffusio n 20 ! 21 X(t) oc t ' 5 . On the other hand, a 
crossover from conventional to anomalous diffusion shows up in the data of system II as the 
points of constant concentration (for which we measure the propagation of density profiles) 
change to higher values. The motion of the point of lowest concentration in system II, 
-^o.oij goes as t ' 5 , and it defines the depth of the Fickian tail ahead of the moving front. A 
slowing-down of the front velocity with time is due to the influence of this Fickian precursor 
as argued by Taylor- 18 . 

Results obtained by a Fourier analysis of the density profiles are shown in Fig. El Starting 
from Fick's equation 

^C(x,t)=D± i C(x,t) i (9) 

where C(x, t) is the density profile in x-direction and D is the collective diffusion coefficient, 
the concentration C(x,t) is represented as a Fourier series (L = L x ) 

C(x,t) = f^C p (t) cos (^p) (10) 



p=0 



C p (t) = y [ L/2 dxC(x,t) COS (^P-) P >0 (11) 
L J-L/2 \ L J 

i r L ' 2 

C (t) = - \ dxC(x,t). (12) 

L J-L/2 

This representation has the advantage that both the periodic boundary conditions and the 
finite system size are automatically taken into account. Then the equation is solved with 
respect to the Fourier coefficients C p , resulting in C p {t) = C p (0) exp(— D(2np/ L) 2 t). It can 
be shown that C p for even p is always as a consequence of our rectangular-shaped initial 
condition; hence we consider here only the odd coefficients. The normalized Fourier coeffi- 
cients £ ?! plotted vs scaled time (2np/L) 2 t should collapse on a master curve (exponential 
decay) if the process was a normal diffusion. The results in Fig. EH follow exactly this be- 
havior and the collective diffusion coefficient estimated from the decay is D pa 0.054. Note 
that using the argument t/L 2 implies the Fickian scaling X ~ t 1 ! 2 . 

In contrast, the behavior of system II (Fig. Eb) reveals a strong deviation from the 
above scenario due to the moving front: The data clearly do not collapse on a single curve 
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(indicating a violation of Fickian scaling), and overshoot into negative values. In the inset, 
a linear scaling X ~ t is assumed. This scaling is better although not perfect. Furthermore, 
it is seen that the reduction of the L exponent "overcorrects" , i. e. the best collapse would 
be observed for an exponent somewhere between 1 and 2. Hence, the data seem to indicate 
that the system is in a crossover region somewhere between the limiting Cases I and II. 

In order to understand what one should expect for the behavior of the Fourier coefficients 
in the limiting Case II, we calculated them for a simple model where a sharp moving front is 
the only feature of the system. We assume that two sharp fronts approach each other with 
constant velocities ±v whereby the profile retains a rectangular shape. One then obtains for 
the odd coefficients C p (t) = C p (0) cos(2npvt / L) . Of course, this result exhibits linear scaling 
X ~ t. Although the corresponding data collapse is not perfect, this result qualitatively 
explains the overshoot into the negative regime, which, according to the analytical result, 
should occur for all modes p > 3. One should keep in mind that the front in system II 
is not a vertical line but has a more complex shape (where different parts of the profile 
follow different laws of motion) and the velocity of the front is not constant throughout the 
simulation but decreases with time. For these reasons, the coefficients show a more complex 
behavior. 

We also tried another type of simulation where the velocity of the front is kept constant 
during the simulation by artificially removing the tails of the profiles: immediately after a 
slow particle is dissolved and is surrounded only by fast particles (the detection of such an 
event being part of the DPD procedure anyway), we remove it from the system and put a 
fast particle in its place. This model attempts to mimic the process of matrix dissolution 
where the slow particles precipitate after they split off from the surface. The set of friction 
coefficients is the same as in system II. In this case the density profile of slow particles looks 
like a step-function with height p and erfc-like boundaries. The front moves linearly with 
time. The Fourier coefficients plotted vs (2ixp/ L)t (data not shown) are then cosine waves 
with same periods (the whole profile shrinks with constant velocity) but different amplitudes 
(due to the erfc-like shape). 

The regime of anomalous diffusion in system II is reached due to the difference in mo- 
mentum relaxation of slow particles in the matrix and in the solvent. In order to check the 
influence of the dynamics of fast particles on the process of dissolution we also use another 
set of friction coefficients: ( ss = 1000, ( s f = 1 and Cff = 0.001. The results (data not 
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shown) are similar to those obtained in system II. This is expected, since even in the limit 
of vanishing friction the interactions in the dense repulsive Lennard- Jones fluid provide a 
particle friction of roughly ( = 20 22 . Whenever Qf is small compared to this value, the 
dynamics of the fast particles is dominated by the conservative forces. 

Eventually, in Fig. 0] we compare our data for the Fourier coefficients C p (t)/C p (0) in 
the case of non-Fickian diffusion to data obtained by using a five times smaller step of 
integration St. Clearly, the results are consistent with one another within the limits of 
statistical accuracy (the curves with the smaller St value were obtained from a single run). 
This demonstrates that the chosen parameters of the MD simulation do not introduce any 
distortion of the physical picture. 

V. CONCLUSIONS 

We have constructed a very simple particle-based computer simulation model which is 
able to exhibit substantial deviations from Fickian diffusion, and shows a clear crossover 
towards Case II behavior (i. e. a constant-velocity front with essentially flat concentration 
profiles ahead and behind). The model disregards all molecular details but incorporates 
the two essential features which we consider as necessary for reaching the Case II limit: A 
substantial discrepancy in the mobility of the two (pure) species, and a "plasticizing" effect, 
i. e. a substantial increase in the mobility of the slow species as soon as its molecules are 
surrounded by the fast species. Our simulation results further corroborate the view that no 
additional molecular ingredients are necessary. We have no definitive explanation why the 
simulations by Tsige and Gres t 20 ^ 21 were unable to reach the Case II limit, but speculate 
that the amount of plasticization might have been insufficient. 

What we find remarkable about our results is not so much the observed physics as such 
- the basic mechanism is essentially rather simple, and the Case II limit should be expected 
as soon as the dynamic asymmetry, combined with plasticization, is strong enough - but 
rather the computational ease with which our model is able to obtain the data, which are 
rather close to the desired limit. The model aims directly at controlling dynamic properties, 
without modifying the equilibrium statics. Furthermore, since we have three parameters 
(Css, Cff an d Csf) & t our disposal, plasticization is particularly straightforward to implement 
by choosing a small value for ( s f. An additional bonus is that hydrodynamics (momentum 
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conservation) is taken faithfully into account. We believe that our approach via DPD is 
particularly promising for simulations of Case II systems. For example, it should now be 
possible to combine our approach with a polymer simulation. There the DPD part would 
model the "glassy" (slow) behavior of the matrix, plus plasticization, while the modeling in 
terms of polymer chains would account for the connectivity and entanglements. Thus in- 
teresting questions of molecular stretching, reorientation, disentanglement kinetics, stresses, 
etc., which are difficult to treat in terms of macroscopic models, could be attacked. 
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Figure captions 

1. Density profiles in x-direction of the slow particles (full symbols) and of the fast 
particles (open symbols) for 4 different times. Part (a): system I, part (b): system II. 

2. (a) Motion of the front position Xi n fi of fast particles for system II. The dashed line 
denotes the dependence characteristic for Case II diffusion X(t) oc t. The data are 
shifted so that their initial position is in the middle of the box X(t = 0.0) = 0.0. A log- 
log plot of X in fi of system II is shown in the inset of part (a) where the slope (dashed 
line) is 0.91 ± 0.01. Three positions of constant density p\ = 0.1 (Xo.i), P2 = 0.05 
(A0.05) and P3 = 0.01 (X .oi) are plotted in (b) against square root of time for the 
slow particles (circles) and for the fast particles (squares). Open symbols denote the 
data of system I, full symbols those of system II. The results are shifted as in part 
(a). Dashed lines in (b) mark the dependence of X(t) typical for Fickian diffusion, 
X(t) oc t - 5 . 

3. The first 9 Fourier coefficients C p (p = 1, 3, 5 ... 17) of the density profiles vs (2Tip/L) 2 t. 
Part (a): system I, part (b): system II. Inset of part (b): The first 9 Fourier coefficients 
C p (p = 1, 3, 5 . . . 17) of the density profiles vs (2np/L)t. 

4. The same as in Fig. Efc> whereby symbols denote data obtained with 5t = 5 x 10~ 5 
while full lines are obtained with 5t = 10~ 5 . 
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